function m = spinMatrix( obj )
% output m is a Structure Array with fields {Imat,Smat,Fmat,IS,F2,mu}
% Imat:1 X 3 cell array giving the nuclear spin matrix in in the direct product space {|I,M_I>|J,M_J>}
%  Smat:1 X 3 cell array giving the nuclear spin matrix in in the direct product space {|I,M_I>|J,M_J>}
% IS: the matrix \vec{I} \dot \vec{J}   
    I2 = obj.I*(obj.I+1); 
    S2 = obj.J.*(obj.J+1);
    
    muI = obj.parameters.mu_I;                 % nuclear magnetic moment in terms of [\mu_{N}]
    LgS = {obj.parameters.LgS, ...              % Lande g-factor for different J
           obj.parameters.LgJ1, ...
           obj.parameters.LgJ2};
    
    Imat=Atom.Spin(obj.I);                     % nuclear spin matrix in the basis {|I,M_I>}
    Smat={Atom.Spin(obj.J(1)), ...             % electron spin matrix in the basis {|J,M_J>}  (1 X 3 cell array)  
          Atom.Spin(obj.J(2)), ...
          Atom.Spin(obj.J(3))};
    
    for i=1:3                                   % opeartors in the direct product space {|I,M_I>|J,M_J>}
        m.Imat{i} = kroneye(Imat, obj.gJ(i));   % Matrix: I X eye(2J+1)  J=1/2,1/2,3/2
        m.Smat{i} = eyekron(obj.gI, Smat{i});   % Matrix: (2I+1) X matrix(J)  J=1/2,1/2,3/2
        m.Fmat{i} = m.Smat{i} + m.Imat{i};      % F=I+J   in direct product space
        m.IS{i} = matdot(m.Imat{i}, m.Smat{i}); % i=1,2 corresponds to J=1/2,i=3 corresponds to J=3/2
        m.F2{i} = (I2+S2(i))*eye(obj.dim(i)) + 2*m.IS{i};
        m.mu{i} = -LgS{i}*muB * m.Smat{i}  + muI*muN/(obj.I+eps) * m.Imat{i}; % muN:nuclear magneton in J/T  referencr:Happer page20-(2.11)
                                                                              % muB:electron magneton in J/T
    end

end

